Critical phenomena in colloid— polymer mixtures: interfacial tension, order parameter, 

susceptibility and coexistence diameter 
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The critical behavior of a model colloid-polymer mixture, the so-called AO model, is studied using 
computer simulations and finite size scaling techniques. Investigated are the interfacial tension, the 
order parameter, the susceptibility and the coexistence diameter. Our results clearly show that the 
interfacial tension vanishes at the critical point with exponent 1v ~ 1.26. This is in good agreement 
with the 3D Ising exponent. Also calculated are critical amplitude ratios, which are shown to be 
compatible with the corresponding 3D Ising values. We additionally identify a number of subtleties 
that are encountered when finite size scaling is applied to the AO model. In particular, we find 
that the finite size extrapolation of the interfacial tension is most consistent when logarithmic size 
dependences are ignored. This finding is in agreement with the work of Berg et al. [Phys. Rev. B, 
47, 497 (1993)]. 

PACS numbers: 61.20. Ja,64.75.+g 



I. INTRODUCTION 

By adding non-adsorbing polymer to a colloidal sus- 
pension, phase separation may be induced. This leads to 
the formation of two coexisting fluid phases, separated 
by an interface 0, ■ The phases are characterized by 
their colloid density, which is high in one phase, and low 
in the other. In order to make the analogy to the fluid- 
vapor transition in atomic liquids, the colloid rich phase 
is usually called the colloidal liquid, and the colloid poor 
phase the colloidal vapor. Obviously, the density of the 
polymers is exactly the opposite: high in the colloidal 
vapor and low in the colloidal liquid. 

In the vicinity of the critical point, a number of impor- 
tant physical quantities are described by simple power 
laws of the form At B . Here, t is some measure of dis- 
tance from the critical point, A is called the critical am- 
plitude, and B the critical exponent To describe 
the critical phase behavior, one thus needs to determine 
the location of the critical point, the critical exponents 
and the critical amplitudes. These quantities can for in- 
stance be obtained in computer simulations, provided fi- 
nite size scaling methods are used. Finite size scaling is 
required because the correlation length diverges at the 
critical point. Since the accessible system size in a sim- 
ulation is finite, the true critical behavior is obscured 
as soon as the correlation length exceeds the size of the 
simulation box 0, Q . Finite size scaling offers a way to 
properly extrapolate the data obtained in simulations to 
the thermodynamic (=infinite system) limit 

In this work, we study the critical behavior of the 
colloid- polymer model introduced by Asakura and Oo- 
sawa [10 (the so-called AO model). In previous simu- 
lations, we have determined the critical point of the AO 
model (for one choice of colloid-to-polymer size ratio) 
and we have provided evidence that this system belongs 
to the 3D Ising universality class [a Q. However, we 
have not yet studied in detail whether we can recover 
the critical exponents, nor have we determined the crit- 



ical amplitudes. The latter quantities are of interest be- 
cause certain critical amplitude ratios are predicted to 
be universal. In this work we combine computer simula- 
tions and finite size scaling to address these issues. The 
quantities that we consider are the order parameter, sus- 
ceptibility, coexistence diameter, and interfacial tension, 
whereby the critical point is approached along different 
paths: from the one-phase region, the two-phase region, 
and along the coexistence line. 

The AO interfacial tension is of particular importance, 
because it gives an indication of the strength of capillary 
waves. This is an issue in (mean-field) density functional 
theories of the AO model. Following Ref.0, the strength 
of the capillary waves is estimated by u = feT '/ '(Aira^ 2 ), 
with a the interfacial tension, £ the correlation length 
in the two-phase region, T the temperature, and fes 
the Boltzmann constant. For 3D Ising critical behav- 
ior, hyperscaling implies that w is constant in the critical 
regime. Note that this constant is universal and given by 
uj ~ 0.8 [ll|. I n contrast, for mean-field critical behavior, 
one would observe a decay of the form uj oc i -1 / 2 . The 
mean-field behavior and the 3D Ising behavior of the cap- 
illary strength are thus profoundly different. Therefore, 
it is important to establish the universality class of the 
AO model. For this purpose, an analysis of the interfacial 
tension is particularly suitable. To determine the criti- 
cal behavior of the interfacial tension, finite size scaling 
methods can be used that do not require prior knowl- 
edge of the universality class. This enables a direct mea- 
surement of the critical exponent, and the corresponding 
critical amplitude. For the AO model, we obtain for the 
interfacial tension a critical exponent 2v = 1.26, which is 
in excellent agreement with the 3D Ising value. 

The critical amplitudes are used to test the universal- 
ity of a number of critical amplitude ratios. Following 
Stauffer 0, 0] , universality implies that only two am- 
plitudes are required to determine the remaining ampli- 
tudes. This allows us to relate the critical amplitudes 
of the AO model obtained in this work, to independent 
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estimates obtained in experimental, theoretical, and sim- 
ulational studies of completely different systems 0, ^| . 
We observe reasonable agreement, but emphasize that 
the error bars, in both our estimates and those in the 
literature, are quite substantial. 

Note that an investigation of the AO critical behavior 
is far more complex than an equivalent study of the 3D 
Ising lattice model would be. For instance, the binodal 
of the AO model is asymmetric, and this gives rise to 
an additional critical power law for the coexistence di- 
ameter. Moreover, the AO model belongs to the class 
of asymmetric binary mixtures, which are generally dif- 
ficult to simulate. In case of the AO model, the accu- 
racy required to apply finite size scaling became avail- 
able only after the recent introduction of a grand canoni- 
cal Monte Carlo cluster move |9j and successive umbrella 
sampling 0] . Note also that the application of finite size 
scaling to asymmetric mixtures [l5| is far less common 
compared to that of symmetric mixtures Q . 

The outline of this paper is as follows. First, the AO 
model is introduced. We then move on to describe the 
simulation techniques used by us. Next, we explain how 
the order parameter, susceptibility, interfacial tension 
and coexistence diameter are extracted from the sim- 
ulation data. The extrapolation of these quantities to 
the infinite system is discussed in section fvTI We then 
present our results and end with a summary in the last 
section. 



II. THE AO MODEL 



The AO model was proposed in 1954 Q and later 
independently by Vrij |7J as a simple description for 
colloid-polymer mixtures. In this model, colloids and 
polymers are treated as spheres with respective radii R c 
and R p . Hard sphere interactions are assumed between 
colloid-colloid (cc) and colloid-polymer (cp) pairs, while 
polymer-polymer (pp) pairs can interpenetrate freely. 
This yields the following pair potentials: 



u cc {r) 

u cp {r) 
Upp(r) 



oo for r < 2R C 

otherwise, 

oo for r < R c - 

otherwise, 



(1) 



= 0, 



with r the distance between two particles. 

Since all allowed AO configurations have zero poten- 
tial energy, temperature plays a trivial role, and the 
phase behavior is set by the colloid to polymer size ra- 
tio q = R p /R c and the fugacities {z c ,z p } of colloids 
and polymers, respectively. The fugacity z a is related 
to the chemical potential [i a via z a — exp(/3/i Q ), with 
a € {c,p}. 

In this work we consider a size ratio q = 0.8 and put 
Rc = 1 to set the length scale. The colloid packing frac- 
tion is defined by rj c = (Air /^)R^N C JV , and the polymer 



packing fraction by rj p = (4tt/3)RIN p /V . Here, iV c (N p ) 
denotes the number of colloids (polymers) inside the sim- 
ulation cell and V the volume of the simulation cell. Fol- 
lowing convention, we use the quantity rf p = z p (47r/3)i?p 
to express the polymer fugacity, rather than z p itself. In 
the literature, rf v is known as the polymer reservoir pack- 
ing fraction. It should, however, not be confused with the 
actual polymer packing fraction in the system rj p . 

The AO model phase separates into a colloidal va- 
por and colloidal liquid , p rovided q and rf p are high 
enough [S [H E I3jpl For q = 0.8, the critical 
point was located at (aQ: 

yyp C1 . = 0.766 ± 0.002, r/ PiCr = 0.3562 ± 0.0006, (2) 
77 CjCr = 0.1340 ± 0.0006, A* CiCr = 3.063 ± 0.003, 

with ^ CjCr the critical value of the coexistence chemical 
potential of the colloids and r] aiCr the critical value of 
r\ a , with a € { p ,p,c}. The above estimates were ob- 
tained using the cumulant intersection method |2l| . and 
by considering field mixing effects p2l l23l | . 



III. SIMULATION METHOD 

We simulate the AO model in the grand canonical en- 
semble. In this ensemble, the fugacities {z Cl z p } and the 
volume V are fixed, while the number of particles inside 
V fluctuates. The simulations are performed in cubic 
boxes with edge length L and using periodic boundary 
conditions. To simulate the AO model efficiently, we use 
a recently developed cluster move 0, . 



A. Phase coexistence 

During the simulation, we measure the probability 
P(rj c ) of observing a certain colloid packing fraction rj c . 
At phase coexistence, the distribution P(r] c ) becomes bi- 
modal, with two peaks of equal area for the colloidal 
vapor and liquid phase. A natural cut-off to separate the 
vapor from the liquid phase is provided by the average 
colloid packing fraction: 



(Vc) = / VcP{Vc)drj c , (3) 
Jo 

where we assume P(rj c ) has been normalized to unity: 

/•oo 

/ P( Vc )A Vc = 1. (4) 
Jo 

The equal area rule simply implies that: 



P(r] c )dr) c 



P{Vc)dr) c 



(5) 



The above equation provides an accurate numerical mea- 
sure to determine phase coexistence [l5j. 
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B. Successive umbrella sampling 

Close to the critical point, the simulation moves back 
and forth easily between the vapor and liquid phases. 
Away from the critical point, at higher polymer fugacity, 
the free energy barrier between the two phases increases. 
In that case, transitions from one phase to the other 
phase become more and more unlikely, and the simula- 
tion will spend most time in only one of the two phases. 
A crucial ingredient in our simulation is therefore the 
use of a biased sampling technique called successive um- 
brella sampling. This technique was recently developed 
by Virnau and Miiller 0] and its purpose is to enable 
sampling in regions where P(rj c ), due to the free energy 
barrier separating the phases, is very low. 



C. Histogram reweighting 

A final ingredient in our simulation is the use of his- 
togram reweighting '2§\. It is based on the observation 
that the probability P(f] c ) measured at one set of model 
parameters (in this case z c and z p ) can be used to esti- 
mate P(f] c ) at different values of these parameters. Ob- 
viously, the gain in computational efficiency is enormous 
because in the ideal case P(r] c ) need only be measured 
once. 

In this work, histogram reweighting is used to locate 
the coexistence fugacity of the colloids. Phase coexis- 
tence is only obtained if the colloid fugacity is chosen 
just right. This value is in general not known at the start 
of the simulation. However, once P(i] c \z c , z p ) has been 
measured at colloid fugacity z c and polymer fugacity z p , 
histogram reweighting can be used to obtain P{rj c \z' c , z^) 
at any other colloid fugacity z' c by using the equation [5( : 

In P(N c \z' c ,z p )= In P(N c \z c ,z p )+(ln^J N c . (6) 

Note that in the above equation rj c has been replaced by 
the number of colloids N c . In the simulations, we thus 
set the colloid fugacity to unity and apply successive um- 
brella sampling to obtain the corresponding probability 
distribution. We then use Eq.© to extrapolate P{rj c ) 
to that colloid fugacity at which the equal area rule of 
Eq.JSJl is obeyed. In practice, it is straightforward to 
write an automated numerical procedure to achieve this. 
We emphasize here that extrapolations in z c are exact, 
in the sense that no statistical or systematic errors are 
introduced (the term "extrapolation" might suggest oth- 
erwise) . 

Within successive umbrella sampling, states (or win- 
dows) are sampled one after the other. In the first win- 
dow, the number of colloids iV c is allowed to fluctuate 
between and 1, in the second window, N c is allowed to 
fluctuate between 1 and 2, and so on. No restriction is 
put on the number of polymers though, so in each window 
N p will fluctuate freely around some average equilibrium 



value (the distribution in N p is to a good approximation 
Poisson like). The crucial point is that by using succes- 
sive umbrella sampling, data over the entire range from 
N c = up to some maximum is obtained (the maximum 
should be chosen well beyond the liquid peak). There- 
fore, extrapolations in z c (which essentially "emphasize" 
the data of some windows with respect to others) can be 
performed without loss of accuracy. 

Histogram reweighting is also used to extrapolate 
P(j] c \z c , z p ) to different polymer fugacities z' p to obtain 
estimates of P(r] c \z c , z' p ). To this end, also the distribu- 
tion of the number of polymers must be recorded for each 
window. Since our implementation of successive umbrella 
sampling only puts a bias on the number of colloids, and 
not on the polymers, this extrapolation will introduce an 
error. The accuracy of the estimated distribution deteri- 
orates when the range z' p — z p over which one extrapolates 
becomes larger. Fortunately, since we are primarily in- 
terested in the behavior close to the critical point, the 
range need not be large and the error introduced by his- 
togram reweighting is small 4]. More importantly, the 
error can easily be checked for as will be shown later. 

Note that in terms of histogram reweighting, the AO 
model is extremely convenient. Since the histograms that 
need to be maintained involve integer data only (namely 
numbers of particles) the problem of choosing a bin size 
for example does not occur. 

IV. EXTRACTING OBSERVABLES 

The coexistence distribution P(rj c ) is a powerful quan- 
tity because a number of important physical observables 
can be extracted from it. For instance, the average pack- 
ing fraction of the colloidal vapor can be written as: 

nine) 

rf c =2 TlcPi^dvc, (7) 

Jo 

and a similar expression holds for the average packing 
fraction of the colloidal liquid: 

/>OG 

ri = 2 / r) c P(r, c )dr) c , (8) 
J (no) 

with (rjc) given by Eq.J3J). In these and following equa- 
tions, the normalization condition of Eq.Q is assumed 
(this also explains the origin of the factors of two in the 
above two equations). 

For the AO model, we define in the two-phase region 
half the difference in colloid packing fraction between the 
vapor and liquid phase as order parameter. The order 
parameter is denoted M c : 

1 V fOC 

M c = ^k=J |r ?c -( 7?c )|P(r ?c )d7 ?C) (9) 

and the coexistence diameter (or rectilinear diameter) is 
given by: 
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where the subscripts "c" emphasize that the definitions 
are expressed in terms of the colloid packing fraction, and 
not the polymer packing fraction. 

In a similar way, we identify the "concentration suscep- 
tibility" or compressibility as the variance of the peaks 
in P(77 c ) |2^. In the two-phase region (if p > we 
define the susceptibility of the colloidal vapor as: 



V 




2 c P(Vo)dr, c 



and the susceptibility of the colloidal liquid as: 



xi = v 




1 \'2 



(11) 



(12) 



where the factors of two are again consequence of the 
normalization condition of Eq.(£j}. Note also the pres- 
ence of the volume V in the above definitions. In the 
one-phase region (r7 p < 7?p cr ), the distribution P(ri c ) will 
loose its bimodal structure and become single peaked. In 
this regime the correct definition for the susceptibility 
reads 0,13: 



V 



r) c P(i] c )dr) c - (r) c ) 



(13) 



with (r] c } given by Eq.©. 

The above definitions are easily modified to define the 
corresponding observables for the polymer phases (sim- 
ply replace the subscript "c" by "p"). Since our sim- 
ulations also store the polymer histograms, coexisting 
packing fractions and susceptibilities can be calculated 
for these phases as well. Note that the above definitions 
are additionally attractive from a numerical point of view 
because the integrations over rj c tend to average out the 
statistical fluctuations that may be present in P(ri c ). 

The interfacial tension is extracted from the logarithm 
of the probability distribution: W = In P(ry c ). Since W 
corresponds to the free energy of the system, the height 
of the peaks in W may be identified as the free energy 
barrier separating the colloidal vapor from the colloidal 
liquid [27]]. In Fig. 2] the barrier is marked Fl, where the 
subscript L emphasizes that the data stem from a finite 
simulation box of size L. In practice, Fl is extracted 
from W via: 



F L = W+-W- 



(14) 



where W+ is the average of W in the peaks: W+ = 
(W v +W{)/2, and W- the value of W at the minimum 
between the peaks (the symbols W v and W\ are defined 
in Fig.^). The corresponding interfacial tension for the 
finite system reads [27| : 



cjl = F L /(2L 2 ), 



(15) 



where the factor of two stems from the use of periodic 
boundary conditions which yield the formation of two 
interfaces in the system. 



Q- 

c 




0.00 0.05 0.10 0.15 0.20 0.25 

FIG. 1: Logarithm of the probability distribution W = 
In P (rjc) for an AO model with q = 0.8 and rf p = 0.765 at 
coexistence. The data were obtained in a cubic simulation 
cell with edge L = 21.0 and using periodic boundary condi- 
tions. The peak at low 77 c corresponds to the colloidal vapor, 
the peak at high rj c to the colloidal liquid. W v (W\) denotes 
the maximum value of the colloidal vapor (liquid) peak. W- 
is the minimum value of W between the two peaks. Fl corre- 
sponds to the free energy barrier separating the vapor phase 
from the liquid phase. 



V. CRITICAL BEHAVIOR 

Essential in the study of critical phenomena is some 
measure of distance from the critical point. As measure 
of distance we use in this work the parameter: 



p,cr 



1 



(16) 



^p.cr) 

This 



which is positive in the two-phase region (rf„ > 
and negative in the one-phase region (?y p < rf v cr ) 
implies that in the one-phase region we must use — t. 

When the critical point is approached from the two- 
phase region, M and a vanish precisely, while \ diverges. 
Common symbols have been established to denote the 
critical exponents and critical amplitudes of the associ- 
ated power laws: 



M = Bt 



J3 



(17) 



where hyperscaling has been assumed (valid for systems 
that belong to the 3D Ising universality class). 

When the critical point is approached from the other 
side, namely the one-phase region, M and a remain zero, 
while x diverges: 

* = r+(-t)-r, (18) 
but with a different critical amplitude r + . 
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The critical behavior of the coexistence diameter D in 
the two-phase region is given by [28|, |2J| : 



D-X, 



At l ~ 



(19) 



with X CI the packing fraction at criticality, given by ry C)Cr 
(??p,cr) hi case of the colloid (polymer) coexistence diam- 
eter. In the one-phase region D is not well defined be- 
cause the distinction between a colloidal vapor and liquid 
is then no longer possible. 

For the 3D Ising universality class, the critical expo- 
nents are given by: 



v 



0.324, 
; 0.629, 



7 : 
a 



1.239, 
0.113. 



(20) 



Also cert ain combin ations of the critical amplitudes are 
Of importance in this work are the 



universal HHHE! 
ratios: 



U 2 = T+/T- 



w 2 Rl' 2 



{Ri) 



3/2 



3/2.p_ 



B 2 



B 2 



4.76 ±0.24, 
3 0.13 ±0.04, 

a 0.71 ±0.13, 



(21) 
(22) 

(23) 



with values taken from Ref. ITTl where k^T is used as 
the unit of energy. The ranges represent the spread in 
different estimates obtained from experiment, simulation, 
and theory. Note that the ranges are quite substantial, 
indicating that the determination of critical amplitudes 
is still a challenge. 



VI. FINITE SIZE SCALING 

In a computer simulation, the critical power laws given 
in the previous section cannot be observed directly. As 
explained before, this is related to the correlation length 
which diverges at the critical point and hence cannot be 
captured in a finite system of size L. However, it is pos- 
sible to perform several simulations at different system 
sizes and use the predictions of finite size scaling theory 
to extrapolate to the thermodynamic limit. Reviews of 
the subject are abundant in the literature 0.EH3 In this 
section, wc will therefore be brief and only reproduce the 
equations required for our analysis. 



A. Finite size scaling of M, \ and D 

According to finite size scaling theory, the order pa- 
rameter Ml obtained in a finite system of linear dimen- 
sion L close to the critical point shows a systematic L 
dependence that can be written as |f|: 



with t the distance from the critical point, critical ex- 
ponents {/3, v} and A4° some function independent of 
system size (A4° is called a scaling function). The criti- 
cal exponents appropriate for the AO model are the 3D 
Ising exponents listed in Eq.J3U|). Since A4° is system 
size-independent, this implies that plots of L@/ v Ml ver- 
sus tL x l v should all collapse onto one master curve, pro- 
vided the correct values of the critical polymer fugacity 
and the critical exponents are used. Such scaling plots |4j 
thus give an indication of whether the assumed univer- 
sality class is indeed correct. Moreover, for large tL 1 ' v 
(but still within the critical region of course) such plots 
should approach the power law behavior of the thermo- 
dynamic limit. This is best visualized if a double loga- 
rithmic scale in the scaling plot is used. The data should 
then approach a straight line, with slope [3 and intercept 
equal to the critical amplitude B. 

The scaling behavior of the susceptibility is analo- 
gously given by: 



XL 



(25) 



In this case, the appropriate quantities for the scaling 
plot are L~*>/ u xl and tL 1 / u . On a double logarithmic 
scale, the data should approach a straight line, with 
slope —7 and intercept equal to the critical amplitude 
r~ or r + (depending on whether the critical point is ap- 
proached from the one-phase or the two-phase region). 

By presenting the simulation data in the form of scaling 
plots, measurements of the critical amplitudes become 
possible. The accuracy of the method increases when 
larger system sizes are used. In order for finite size scal- 
ing theory to apply, it is important that the simulations 
are performed in the so-called scaling regime. At which 
system size L the scaling regime begins, varies from sys- 
tem to system and is a priori not known. If, however, the 
considered system sizes are too small, it will show in the 
corresponding scaling plots as systematic deviations from 
any data collapse onto a master curve. Alternatively, one 
could consider corrections to finite size scaling theory for 
these smaller systems |30| . 

In principle, a finite size scaling analysis of the coex- 
istence diameter D is also possible. In this case, the 
appropriate scaling plot is (D — X cx )L( 1 ~ a " v expressed 
as a function of tL 1 ^ . In practice, however, it is diffi- 
cult to distinguish in simulation data the term t 1_a in 
Eq. I|19(l from the next order term in the series expansion 
(which would be linear in t) because the critical exponent 
a is rather small. The accuracy of such scaling plots is 
therefore usually rather poor. 



B. Finite size scaling of a 

The interfacial tension in the thermodynamic limit <7oo 
(in d dimensions) is related to the free energy barrier of 
the finite system Fl via 27]: 



Mr 



L-^ v M°(tL^ v ), 



(24) 



exp(F L ) = AL x exp(2L d - L a 00 ) 



(26) 
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where trivial factors of fc^T have been dropped. Taking 
the logarithm on both sides, the above equation can be 
written as: 



;1nL 



In A 



2L d " 1 2L d ~ 1 ' 



(27) 



with a i, the interfacial tension of the finite system given 
by Eq.JEJ and constants that are generally not 

known. While it is not possible to measure Coo directly, 
it is possible to measure the interfacial tension of the 
finite system <jl for several system sizes L, and then use 
the above equation to extrapolate to the thermodynamic 
limit. 

One attractive feature of Ea. (|27[) is that it does not 
depend on the critical exponent v and can therefore be 
applied without prior knowledge of the universality class 
of the system. In other words, it can be used to measure 
v as is demonstrated in Ref. for the Lennard-Jones 
fluid. For the AO model, the universality class is already 
known 0, Q . In this case, Eq. I|27(l still provides a power- 
ful consistency check: if the universality class of the AO 
model is indeed 3D Ising, it should be possible to extract 
the exponent v from the simulation data. 

The issue here is how to perform the extrapolations in 
practice |32t I33L 13 ll | . Ideally, the simulation data for the 
different system sizes should be extrapolated to the ther- 
modynamic limit using two-parameter fits in the vari- 
ables ln(L)/L d_1 and \jL d ~ x . This approach, however, 
may have numerical problems associated with it. Typi- 
cally, only data over a relatively small range of different 
system sizes is available (this is certainly the case for a 
non-trivial mixture like the AO model). It will be dif- 
ficult to separate the \a{L)/L d ~ 1 term accurately from 
the \jL d ~ x term over such a small range. Therefore, this 
extrapolation scheme will likely lead to poor precision. 

Alternatively, one can argue that for small L it may 
occur that a;lnL| < |ln^4|, in which case an extrapola- 
tion in the single variable is most appropriate, 
whereas for large L the si ngle variable \a(L)/L d ~ 1 is the 
better choice (27J . In Ref. |3l], for example, the extrapo- 
lations are performed in the single variable \n.(L)/L d ~ l . 
Since it is a priori not clear which of the above extrapo- 
lation methods is the most accurate, all are investigated 
in this work. 



VII. RESULTS 

The following results stem from simulations of the AO 
model with q = 0.8 performed in cubic boxes with edge 
length L and using periodic boundary conditions. The 
dimensionality of the simulations is d = 3. In order to 
apply finite size scaling, the following system sizes are 
considered: L\ — 15.5, Li = 16.7, L 3 = 17.7 and L4 = 
21.0. For each system size, the coexistence probability 
P(tj c ) is measured accurately at one value of ?7 p chosen in 
the vicinity of the critical point. Histogram reweighting 
is used to extrapolate P{j] c ) to other values of rf p . We also 




ln(x = tl_' 



FIG. 2: 3D Ising scaling plot for the order parameter of the 
colloids M c given by Eq.® for the AO model with q = 0.8 in 
the two-phase region. The choice of axis is explained in the 
text. The critical amplitude is extracted by means of a fit to 
the tail of the data, see also Table □ 



performed a number of shorter simulations to explicitly 
measure P{rj c ) at different values of rf p . The results of 
these simulations were used to check the consistency and 
accuracy of the extrapolated distributions. The quality 
of our data is such that extrapolations over the range 
rf p w 0.72 to ?7p ss 0.83 can be carried out reliably. 



A. Order parameter and binodal 

The critical behavior of the order parameter of the col- 
loids M c is analyzed in the scaling plot of Fig. [5] The 
collapse of the simulation data from the different system 
sizes onto one master curve is clearly visible. The scaling 
plot for the order parameter of the polymers M p is shown 
in Fig. |21 which again demonstrates the collapse onto one 
master curve. In these figures, the critical polymer fu- 
gacity ?7p cr was used as a free parameter and tuned until 
the best collapse occurred. By performing a linear least 
squares fit to the tails of the master curves, the critical 
amplitudes can be obtained. The corresponding critical 
power laws are given in Table U which also lists the value 
of ?7 pjCr that was used in the scaling plots. Naturally, 
this value of rf p cl . should agree with the previous esti- 
mate listed in Eq.©. The critical amplitudes obtained 
from the fits are sensitive to the range over which the fit 
is performed: the variation is used as a measure for the 
error in Table [I] 

The critical behavior of the coexistence diameter of 
the colloids D c is presented in Fig. 21 In this case, the 



best collapse of the data occurs at ryl, 



0.771 which 



7 





FIG. 4: Scaling plot of the coexistence diameter of the col- 
loids D c . The value of X CI — r] CiCr used in this plot was taken 
from Eq.(|2J. Allowing variations in X cr did not improve the 
collapse of the data. 



is not in agreement with Eq.J2J). Some discrepancy is 
to be expected though because the singularity in the co- 
existence diameter is very weak, which makes it hard 
to discern it from simulation data. The behavior of the 
data in the tails, however, seems rather well described by 
the exponent 1 — a. The corresponding scaling plot for 
the coexistence diameter of the polymers D p is shown in 
Fig. [SJ In this case, the curvature of the data for small t 
seems in error. Therefore, we conclude that our simula- 



FIG. 6: Phase diagram of the AO model in reservoir repre- 
sentation. The solid curve (L = oo) shows the binodal in the 
thermodynamic limit obtained using the power laws for M c 
and D c listed in Table Q] (here rjL cr = 0.766 was used). The 
dashed/open curves are raw simulation data for various finite 
system sizes L as indicated. The bar marks the location of 
the critical point given by Eq.@. 



tion data is not accurate enough to reliably extract the 
critical behavior of the coexistence diameter. The power 
laws for D c and D p , given for completeness in Table [I] 
must therefore be treated with some care. 

By combining the expressions for M c and D c in Ta- 
ble [U the colloid packing fractions of the vapor and liq- 
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FIG. 7: Phase diagram of the AO model in system repre- 
sentation. The solid curve (L = oo) shows the binodal in the 
thermodynamic limit obtained using finite size scaling (here 
^p,cr — 0.766 was used). The dashed/open curves are raw 
simulation data for various finite system sizes L as indicated. 
The square marks the location of the critical point given by 
Eq.@. 



uid phase (rjl and r)\) can be written as functions of rf v . 
These expressions, which are valid close to the critical 
point in the thermodynamic limit, describe the binodal 
of the AO model in reservoir representation. The result- 
ing binodal is shown in Fig. together with the raw 
simulation data from the various system sizes. The fig- 
ure clearly shows the familiar finite-size deviation of the 
simulation data close to the critical point . If also the 
critical expressions for M p and D p are used, rf p can be 
eliminated altogether to yield the binodal in the experi- 
mentally more relevant (j] c ,Vp) or system representation. 
The resulting plot is shown in Fig. [7| together with the 
raw simulation data. 



B. Susceptibility 

Next, we consider the critical behavior of the suscep- 
tibility. Fig. [S] shows the scaling plot of the suscepti- 
bility of the colloidal phase Xc in the one-phase region. 
The critical power law extracted from this plot is listed 
in Table |U The scaling plot for the susceptibility of the 
polymer phase Xp m the one-phase region is qualitatively 
similar and not shown. Instead, only the critical power 
law is given in Table [I] 

Measurements of the susceptibility in the two-phase 
region are prone to a number of potential numerical pit- 
falls. Since the susceptibility is a second order moment 
of the distribution P{rj c ), it is generally more sensitive to 
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FIG. 8: Scaling plot of the susceptibility of the colloidal phase 
Xc in the one-phase region. In this case Eg. (1131 is used to 
measure \c- 



statistical errors than first order moments like the order 
parameter. For asymmetric systems like the AO model, 
the additional problem arises that the statistics in the 
liquid peak of P(r) c ) are systematically better than in 
the vapor peak: simply because the liquid peak contains 
more colloids. These problems will of course vanish with 
increasing system size. However, for a non-trivial system 
like the AO model, the system sizes that can be handled 
today are unfortunately still in the regime where these 
subtleties come into play. 

In case of the colloidal vapor and liquid in the two- 
phase region, the relevant susceptibilities are xZ given by 
Ea. (|ll|) and x c given by Ea. (|12fl . In principle, close to 
the critical point, the susceptibility should be the same 
in both phases. In practice, there will be differences due 
to the numerical difficulties outlined above. Two proce- 
dures are now conceivable: (1) use the average of xZ an d 
Xc a s best estimate of the susceptibility, or (2) use only 
xi- Which of these procedures is the best needs to be 
checked with the data available. In our case, the suscep- 
tibility of the colloidal phase in the two-phase region was 
best described using only x c - The resulting scaling plot 
is shown in Fig. Eland the corresponding power law in Ta- 
ble[I] The collapse of the data is clearly visible, but unlike 
Fig. the slope of the data for large t is not quite —7. 
A less satisfactory fit is expected though, because the 
statistics in Fig. |5| are significantly worse compared to 
Fig. [S] since only half the data is used. 

The susceptibility measurements of the polymer phases 
in the two-phase region were most accurate if the average 
Xp 1 = (xp + Xp)/2 was used. The scaling plot of xp 1 
is shown in Fig. I1UI and the corresponding power law is 
listed in Table |U As in Fig. [51 the data collapse is clearly 
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FIG. 9: Scaling plot of the susceptibility of the colloidal liquid 
Xc given by Eg. 1121 in the two-phase region. 
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FIG. 10: Scaling plot of the average susceptibility of the two 
polymer phases Xp m the two phase region. 
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FIG. 11: Susceptibility of the colloids Xc across the phase 
transition. The solid curve shows the susceptibility in the 
thermodynamic limit given by the critical power laws of Ta- 
ble m with 77p, cr = 0.766. The remaining curves show the raw 
simulation data of the various system sizes L. 



C. Interfacial tension 

As mentioned before, the finite size extrapolation of 
the interfacial tension is less straightforward and vari- 
ous methods can be used. Since the number of differ- 
ent system sizes considered by us is rather small, multi- 
parameter fits such as the ones investigated in Ref. |3j| 
are out of the question. Instead, we investigate single 
parameter fits only, which in this case essentially means 
choosing \n(L)/L d ^ 1 or 1/L d ~ 1 as scaling variable. The 
results of both extrapolation procedures for a number 
of different rf p are summarized in Fig. El and Fig. IT31 
The extrapolations produce meaningful estimates up to 
rjp « 0.80. For higher values of 77^, one leaves the criti- 
cal regime and Eg. 1)27(1 breaks down. In this regime, the 
interfacial tension gradually becomes system size inde- 
pendent. 

The fits in Fig. Upland Fig. El seem equally accurate, 
and based on these figures alone we cannot reject one 
extrapolation method in favor of the other. The issue is 
resolved when the interfacial tension Uoo in the thermo- 
dynamic limit itself is considered. In this case, the critical 

2v 



demonstrated, but the slope -7 is not quite reproduced. P ower law ^oo = ero* is valid which implies that plots of 



In Fig. ^2 we plot the susceptibility of the colloids as 
function of rf p in the vicinity of the critical point. Shown 
are the raw simulation data as well as the critical power 
laws of Table [I] Clearly demonstrated is the familiar 
finite-size rounding of the raw simulation data in the 
vicinity of the singularity Q. 



(Too versus t on double logarithmic scales should collapse 
onto straight lines, provided the correct value of rj^ cr in 
t is used. Note that ?7p iCr is the only free parameter: the 
critical exponent v follows automatically from the slope 
of the line. The resulting plots are shown in Fig. 1141 and 
Fig- El i n which ln(L)/L <i_1 and l/L^ -1 were used as 
scaling variable, respectively, and rf v cr in each plot was 
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FIG. 12: Finite size extrapolation of the interfacial tension 
in ln(L)/L d_1 . This extrapolation scheme corresponds to the 
assumption ln^4 = in Eg. 1271 . Shown in the above plot is 
the interfacial tension of the finite system ox as function of 
ln(Z/)/I/ d-1 for various values of rf p as indicated. The straight 
lines are linear least squares fits to the data. The intercept 
of these lines with the ordinate yields an estimate for the 
interfacial tension Ooo in the thermodynamic limit. 



tuned until the best collapse occurred. In these figures, 
measurements of the interfacial tension up to rf p = 0.79 
were used. 

One important observation is that both data sets in 
Fig. ^] and Fig. accurately reproduce the expected 
slope 2v w 1.26 corresponding to 3D Ising behavior. 
However, if \n.(L)/L 1 is used as scaling variable, the 
best collapse is obtained at rf cr = 0.7696, which is not 
in agreement with the previous estimate of Eq.Q. On 
the other hand, if l/i d_1 is used, the best collapse is ob- 
served at ?7p cr = 0.7661, which is in excellent agreement 
with Eq.©. Therefore, we conclude that 1/L d_1 is the 
appropriate scaling variable for our problem. The cor- 
responding power law obtained from Fig. 1151 is listed in 
Table H] The effect of using the incorrect scaling variable 
seems to be that the critical point is not correctly esti- 
mated, while the critical exponent is less affected. It thus 
seems wise practice to always compare the critical point 
obtained from finite size extrapolations of the interfacial 
tension to some other independent estimate (this esti- 
mate could for instance be obtained using the cumulant 
intersection method). 

It is worth noting that Berg et al. in Ref. also con- 
clude that the extrapolation of the interfacial tension is 
most consistent if x = is assumed in Ea. (|27|l . In ret- 
rospect, it is not surprising that l/U 1 ^ 1 is the appro- 
priate scaling variable in our case. Close to the critical 
point, the distribution P{rj c ) scales with the system size 



FIG. 13: The analogue of Fig.[T21but using 1/L d ~ 1 as scaling 
variable instead. This extrapolation scheme corresponds to 
the assumption x — in Eg. 1271 . 



as unmni: 

P L (ri c ) = b L^r (b L^ri c ), (28) 

where Pl{Vc) is the distribution P(r) c ) measured in the 
finite system of size L, bo some non-universal constant, 
and V° a function independent of system size. According 
to Ea. l|14|l . the free energy barrier Fl is given by the 
peak-to-valley height in the logarithm of Pl(t)c)- The 
above scaling property thus implies that Fl becomes L 
independent close to the critical point. As a result, the 



TABLE I: Summary of the critical behavior of the AO model 
with q — 0.8. Shown are the critical power laws in the one- 
phase and two-phase regions of various physical quantities 
obtained in scaling plots. Listed in the last column is the 
value of 7/p, C r at which the best collapse in the corresponding 
scaling plot was observed (provided here to give an indication 
of the consistency of our results). The symbols M, D, \ an< l 
a are defined in section ITV1 while the critical exponents are 
listed in Eo.CT. 
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FIG. 14: Interfacial tension in the thermodynamic limit 
(Too as function of t on a double logarithmic scale, where 
ln(L)/L d_1 was used as scaling variable. The best collapse 
onto a straight line is observed at rf p ^ cr = 0.7696, which is not 
in agreement with Eq.Q. The slope of the line yields the 
critical exponent of the interfacial tension, for which we find 
2u = 1.25 ±0.01. 
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FIG. 15: The analogue of Fig. EU in which l/L d ~ x is 
used as scaling variable instead. The best collapse occurs 
at ?7p,cr = 0.7661, which is in excellent agreement with Eq.(|5J. 
For the critical exponent we obtain 2v = 1.26 ± 0.01 and for 
the critical amplitude ctq = 0.26 ± 0.02. 



FIG. 16: The analogue of Fig. 1141 using two parameter fits 
in ln(L)/I/ d_1 and 1/L d ~ 1 when extrapolating the interfacial 
tension to the thermodynamic limit. The best collapse occurs 
at ^p.cr — 0.7689. For the critical exponent we obtain 2v = 
1.30 ± 0.04 and for the critical amplitude oo = 0.33 ± 0.03. 



interfacial tension should scale with 1/L d_1 . 

For completeness, we also show in Fig. ^| the behav- 
ior of (Too as function of t when two-parameter fits in 
both ln(L)/L d_1 and \/L d ' 1 are used. Since only four 
different system sizes are considered, the statistical qual- 
ity of the extrapolations will likely deteriorate, and this 
clearly shows in Fig.^| Still, even then, it is possible to 
extract the critical exponent and the critical point with 
reasonable accuracy. Note that the deviations in 77^ cr ob- 
tained using the various extrapolation schemes are below 
the one percent level. In practice, the choice of extrap- 
olation scheme thus only becomes important when high 
resolution data are available. 



D. Critical amplitude ratios 

We now turn to the critical amplitude ratios that can 
be extracted from Table [I] We first calculate U2 given 
by Eq. I|21|l . If we consider the colloidal phases we obtain 
U2 = 3.9 ± 0.9. The corresponding ratio for the polymers 
is found to be U2 = 4.1 ± 1.9. Within the error bars, 
these estimates are compatible with Ea. H21|) . The quan- 
tity w 2 R 3 J 2 is found to be 0.10 ± 0.04 for the colloidal 
phases, and 0.08 ± 0.04 for the polymer phases. This 
is again compatible with Eci. (|22)l . although one must be 
aware of the large error bars in our estimates. Finally, 
the quantity {R+f/ 2 /Q c is found to be 0.40 ± 0.16 for 
the colloids, and 0.35 ± 0.07 for the polymers. In this 
case, we systematically underestimate the values listed 
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in Ea. l|23[) . Note, however, that the discrepancy with the 
3D Ising values is not too severe. Given the difficulty in 
general of measuring critical amplitudes, even in the case 
of simple lattice models 0] 1 the agreement we obtain is 
already quite remarkable. 

VIII. SUMMARY AND OUTLOOK 

In summary, we have studied the critical behavior of 
the order parameter, interfacial tension, susceptibility 
and coexistence diameter of the AO model with colloid 
to polymer size ratio q = 0.8. An important result is that 
the critical exponent of the interfacial tension equals the 
expected 3D Ising value 2v « 1.26. This critical behav- 
ior is consistent with previous simulations H,i|, and also 
with experimental work (3l} in which 3D Ising critical 
behavior was observed in real colloid-polymer mixtures. 
The critical behavior of the order parameter and the sus- 
ceptibility are also 3D Ising like. Our data for the co- 
existence diameter is not conclusive. This is related to 
the small value of the critical exponent a, which makes it 
difficult to accurately resolve the critical behavior from 
the simulation data. More accurate simulations of larger 
systems are required to resolve this. 

We found that the critical amplitude ratios obtained 
from our simulations are compatible with the 3D Ising 
universality class. This confirms the consistency of our 
data, and is encouraging considering the AO model is 
an asymmetric binary mixture and therefore difficult to 
simulate. We emphasize, however, that our estimates for 
the critical amplitudes should be regarded as consistency 
checks only. Their accuracy cannot compete with that 
obtained in, for example, direct simulations of the 3D 



Ising lattice model. 

We have demonstrated that finite size scaling meth- 
ods are equally applicable to rather complex systems like 
the AO model, and can be used to obtain results with 
meaningful accuracy. Regarding finite size extrapolations 
of the interfacial tension, our analysis is consistent with 
Ref. |U in the sense that the most consistent fits are ob- 
tained by ignoring the ln(L) dependent term in Ea. (|27|l . 
For the AO model, we also observed that by using the 
incorrect scaling variable, the critical point in particular 
is not estimated correctly. 

This work provides additional insight into the critical 
behavior of the AO model. For a complete understand- 
ing, the critical behavior of the correlation length should 
still be investigated. While the corresponding critical ex- 
ponent is of course known, namely —v, the critical am- 
plitude is not. This critical amplitude, however, can- 
not be extracted from the probability distribution P(rj c ). 
The usual approach to study the correlation length is to 
consider the static structure factor instead [26j. This re- 
quires additional simulations which are beyond the scope 
of this work, and will be postponed to a future publica- 
tion. 
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